#include "solarposition.h"

#include "globals.h"
#include "util.h"
#include <cmath>

namespace libsgp4
{

Eci SolarPosition::FindPosition(const DateTime& dt)
{
    const double mjd = dt.ToJ2000();//计算J2000.0历元的儒略日MJD
    const double year = 1900 + mjd / 365.25;//近似年份，用于天文计算的时间基准
    const double T = (mjd + Delta_ET(year) / kSECONDS_PER_DAY) / 36525.0;//从J2000.0历元开始的时间参数，单位为世纪（Julian Centuries）
    const double M = Util::DegreesToRadians(Util::Wrap360(358.47583
                                                          + Util::Wrap360(35999.04975 * T)
                                                          - (0.000150 + 0.0000033 * T) * T * T));//太阳的平近点角 M，描述太阳在轨道上的平均位置
    const double L = Util::DegreesToRadians(Util::Wrap360(279.69668
                                                          + Util::Wrap360(36000.76892 * T)
                                                          + 0.0003025 * T*T));//太阳的平黄经 L，描述太阳在黄道上的平均位置
    const double e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T;//地球的轨道偏心率，描述地球轨道的椭圆形状
    const double C = Util::DegreesToRadians((1.919460
                                             - (0.004789 + 0.000014 * T) * T) * sin(M)
                                            + (0.020094 - 0.000100 * T) * sin(2 * M)
                                            + 0.000293 * sin(3 * M));//太阳的章动修正项，用于修正平近点角到真近点角
    const double O = Util::DegreesToRadians(
        Util::Wrap360(259.18 - 1934.142 * T));//太阳的近地点角距，用于修正黄经
    const double Lsa = Util::WrapTwoPI(L + C
                                       - Util::DegreesToRadians(0.00569 - 0.00479 * sin(O)));//太阳的真黄经 Lsa
    const double nu = Util::WrapTwoPI(M + C);//太阳的真近点角，表示太阳相对于地球在椭圆轨道上的位置。
    double R = 1.0000002 * (1 - e * e) / (1 + e * cos(nu));//太阳到地球的距离 R，单位是天文单位AU
    const double eps = Util::DegreesToRadians(23.452294 - (0.0130125
                                                           + (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * cos(O));//黄赤交角
    R = R * kAU;

    Vector solar_position(R * cos(Lsa),
                          R * sin(Lsa) * cos(eps),
                          R * sin(Lsa) * sin(eps),
                          R);

    return Eci(dt, solar_position);
}

double SolarPosition::Delta_ET(double year) const
{
    // 计算地球时间（ET）与协调世界时（UTC）之间的差值（单位：秒）
    // 公式包含三部分：
    // 1. 常数偏移量 26.465，用于修正基准时间差。
    // 2. 线性变化项 0.747622 * (year - 1950)，表示时间差随年份的线性变化。
    // 3. 周期性修正项 1.886913 * sin(kTWOPI * (year - 1975) / 33)，考虑地球轨道周期性变化的影响。
    // 其中，1950 年是线性变化的参考年份，1975 年是周期性变化的参考年份，周期为 33 年。
    return 26.465 + 0.747622 * (year - 1950) + 1.886913
                                                   * sin(kTWOPI * (year - 1975) / 33);
}

} // namespace libsgp4
